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The principle objectives of NCC3-885 

• Validate, access scalability, and apply two parallel flow solvers for three-dimensional 
Navier-Stokes flows. 

• Develop and validate high-order parallel solver for Direct Numerical Simulations 
(DNS) and Large Eddy Simulation (LES) problems. 

• Investigate and develop high-order Reynolds averaged Navier-Stokes turbulence 
model. 

Validate, Access Scalability, and Apply Two Parallel Flow Solvers For Three-Dimensional 
Navier-Stokes Flows 

Two parallel flow' solvers were validated and used for investigating basic fluid flows. Both 
parallel solvers are versions of an implicit three-dimensional time-accurate Beam-Warming 
Navier-Stokes solver. Both are second-order accuracy in time and space. 

The first parallelization methodology modified a flow solver to decompose a single block grid 
and the solver algorithm using two-dimensional multi-partitioning which evenly distributed the 
work across multiple processors. This approach modified the solver algorithm such that the 
multi-partitioned grid is solved as if it were still a single domain. An example of a single domain 
decomposition for four processors is seen in Figure 1. Inter-grid flow information is 
communicated between the processors using the Message-Passing Interface (MPI) library. The 
scalability of this approach maintained about 50% efficiency for 64 processors on two 
supercomputers, the IBM SP2 and the Silicon Graphics 02K, as shown in Figure 2. The code 
was validated for Couette flow and steady/unsteady flow over a circular cylinder. Two and three- 
dimensional course and fine grid simulations were accomplished for both stationary and rotating 
circular cylinders. A comparison of the flow for the 2-D for the stationary and rotating cylinder 
is shown in Figures 3 and 4. Using cylinder rotation as a form of flow control resulted in a 69% 
(2-D) and 41% (3-D) reduction in the mean drag when compared to the stationary counterpart. 

Results of the above cases were presented in a paper [1] at the AIAA 38th Aerospace Science 
Meeting & Exhibit at Reno, NV in January 2000. Investigations of these flows were expanded 
and incorporated into a journal paper [2], 



Figure 1. Multi-partion scheme for single block grid. 


Figure 2. Scalability of Multi-Partitioning Code. 



Figure 3. Z-vorticity contours 2-D Cylinder at 
Re=3900. 


Figure 4. Z-vorticity contours 2-D Rotating Cylinder 
at Re=3900. 


The second approach parallelized a multi-block overset grid domain decomposition (Chimera) 
scheme using the MPI library for information communication between the distributed grids. The 
new solver was validated for Couette flow and both steady and unsteady flow over a circular 
cylinder. The code was then applied to flow over an ogive-cylinder (Figure 5)and a finned 
missile (Figure 6). Scalability, shown in Figure 7, of the Chimera based approach proved 
superior to the first methodology since the parallel solver was able to maintain 90% efficiency 
for up to 64 processors on two supercomputers, the IBM SP3 and the Silicon Graphics 02K. 
Results from this research were presented at the AIAA 39th Aerospace Sciences Meeting and 
exhibit in Reno, NV in January 2001 [3]. 
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Figure 5. Total Pressure Contours and Streamlines Figure 6. Streamlines and Total Pressure Contours for 
for Ogive-Cylinder. Finned Missile at 20 degrees AOA. 



Figure 7. Speed-up for Parallel Chimera based solver 

The parallel overset grid based solver was subsequently used to study the flow over a 3-D 
pitching NACA 0012 wing section attached to an end-wall being pitched at a constant rate from 
zero incidence to 60 degrees angle of attack. Two-dimensional solutions were generated as a 
baseline for comparison to 3-D results as well as results from a 2-D sixth-order scheme and 
previous numerical and experimental work. The location and angles when the primary, 
secondary, and tertiary leading-edge vortex formed in the 3-D flow were in close agreement with 
the 2-D second-order and sixth-order schemes. The inception of the leading edge vortices for 3- 
D calculations were virtually unaltered from the 2-D findings. 

Three-dimensional effects were found to originate in the leading edge secondary flow structures 
after the development of the three leading edge vortices. Breakdown of the primary dynamic stall 
vortex started at the wall and propagates to the symmetry plane over eleven degrees of the pitch- 
up maneuver. After the onset of three-dimensional effects, the flow features of the 3-D symmetry 
plane and the 2-D solutions departed as the pitch angle increased. 
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Figures 8 and 9 show 2-D flow features at various stages of the dynamic stall process. Figures 10 
and 1 1 display the development of cellular secondary flow structures at the leading edge and 
breakdown of the dynamic stall vortex, respectively. These results were presented at the 31st 
ALA. A Fluid Dynamics Conference and Exhibit at Anaheim, CA in June 2001 [4]. 




Figure 8. Streamlines, Vorticity Contours, and Surface 
Pressure at 19.68 degrees AOA. 


Figure 10. Z-Vorticity Contours for 
32.03 degrees AOA. 



Figure 9. Streamlines, Vorticity Contours, and Surface 
Pressure at 40. 10 degrees A. 


Figure 11. Iso-Vorticily Surface and Z-Vorticity 
Contours for 34.05 degrees AOA. 


Develop and Validate High-Order Parallel Solver For Direct Numerical Simulations (DNS) 
And Large Eddy Simulation (LES) Problems 

The next major effort consisted of incorporating the parallel overset grid based approach into a 
high-order compact-difference solver for application to LES and DNS problems. The new solver 
provides temporal second-order implicit and fourth-order explicit time-integration schemes and 
up-to sixth-order spatial accuracy. Numerical stability is maintained using a tenth-order low-pass 
filter. Validation of the new solver was successfully accomplished using three benchmark cases: 
unsteady Couette flow, and inviscid convecting vortex, and flow over a circular cylinder. 
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Application of the new solver to large eddy simulation of turbulent channel flow and flow past a 
cylinder with a turbulent wake demonstrated the ability of the code to accurately and efficiently 
simulate fine-scale fluid structures. Turbulent channel flow simulations were performed for the 
Reynolds numbers of Rex=180, 395, and 590. For the case of Rex=180, the flow was solved both 
with and without a dynamic sub-grid scale model. The no-model approach is also known as the 
Implicit LES (ILES) methodology. The ILES method was able to match DNS results as well as 
the simulations with the dynamic SGS model at half the computational cost. The channel flow 
results at Rex=395 and 590 also agreed well with DNS solutions. 

The final case investigated to validate the code during this initial phase of development was 
transitional flow over a circular cylinder at ReD=3900. This study started by directly comparing 
the vector and parallel versions of the high-order compact-difference code using identical initial 
conditions, boundary conditions and solver parameters. The vector and parallel solvers resulting 
mean flow quantities compared well. The flow over the cylinder was also solved on a fine 12 
million point grid. The flow characteristics of the fine grid solution matched the experimental 
mean streamwise velocity profile in the near wake of the cylinder which previous coarse grid 
solutions and other numerical studies were unable to properly obtain. 

Figures 12 and 13 demonstrate the ability to resolve fine-scale fluid structures with the new 
parallel high-order solver. These results were presented at the 32nd AIAA Fluid Dynamics 
Conference and Exhibit at St. Louis, MO in June 2002 [5], 



Figure 12. Instantaneous streamwise velocity contours at Figure 13. Instantaneous midspan z-vorticity 

channel inllow an at y+=5.46. contours for a circular cylinder. 


The next application of the parallel high-order compact-differencing solver was to a large eddy 
simulation of a Aerospatiale A-airfoil wing section at near stall conditions and Re=2. 1x106. 
Parameters were investigated including impact of the spanwise and global grid resolution, order 
of spatial discretization, trailing edge geometry, angle-of-attack, and use of a SGS model. 
Resulting from this parametric study was that improved spatial resolution and scheme accuracy 
had the primary influence on properly resolving flow structures and capturing transition. 
Improved spatial resolution and accuracy of the scheme resulted in increased levels of turbulence 
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which decreased the required computational spanwise extent. The high-order discretization 
displayed the ability to capture natural transition even on fairly coarse grids without introducing 
special scheme modifications. All cases computed in this work correctly captured the natural 
transition near the leading edge of the airfoil within two-percent chord length of the experimental 
value. 

Figures 14 and 15 show the difference that the spatial discretization order make in solution 
results. Figure 14 shows a standard second-order discretization resulting in very little small-scale 
flow structure resolution and Figure 15 with significantly better modeling of the flow physics. 
These results were presented at the AIAA 41st Aerospace Science Meeting & Exhibit at Reno, 
NV in January 2003 [6] and at the 33rd AIAA Fluid Dynamics Conference and Exhibit at 
Orlando, FL in June 2003 [7, 8], 



Figure 14. U-velocity contours for 
second-order scheme 



Figure 15. U-velocity contours for 
high-order scheme 


Investigate and Develop High-Order Reynolds Averaged Navier-Stokes Turbulence Model 

The final effort during this cooperative agreement was development of a high-order k-e 
turbulence model. The model was validated for flow over a flat plate and compared with a 
second-order version of the modei for both second- and nigh-order flow solvers. This portion of 
the study demonstrated that a second-order turbulence model dominates the solution even when 
high-order compact differencing is used for the flow equations. A grid resolution study verified 
the expected trend that the solution reached an asymptotic state much faster with a high-order 
turbulence model that with a standard second-order model as shown in Figures 16 and 17. 

The high-order k-e turbulence model was then applied to flow over a wall-mounted hump. 
Simulations for the hump were completed for conditions with and without flow control and were 
compared to previous numerical and experimental results as part of the NASA Langley Research 


6 



Figure 16. Turbulent Flat Plate Velocity Profile Figure 17. Turbulent Flat Plate Velocity Profile 

for Coarse Grid. for Fine Grid. 

Center CFD Validation on Synthetic Jets and Turbulent Separation Control Workshop held in 
Williamsburg, VA in March 2004. Some disagreement with experimental data was seen in the 
current study. This was also seen by all the members that studied this configuration for the 
workshop. Wind tunnel blockage is believed to be responsible in part for the disagreement. 
Comparison of the compact solver with a high-order turbulence model and a second-order 
flow/turbulence solver were also performed. Qualitative agreement was achieved with 
experimental data for both high- and low-order schemes. High-order solutions on a coarse grid 
agreed very well with second-order solutions on a grid four times larger. The reattachment 
location in the numerical simulations was consistently less than 10% further downstream than 
experimental values. Figures 18 and 19 display comparisons of u-velocity contours and 
streamlines between experimental PIV and numerical results. 



Figure 18. Comparison of Experimental PIV (top) and Figure 19. Comparison of Experimental PIV (top) and 
Numerical (bottom) U-velocity Contours. Numerical (bottom) Streamlines. 
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Both the flat plate flow and wall-mounted hump cases demonstrated that significant net 

computational savings are possible when using high-order discretization for the flow equations 

and the turbulence model. 
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Dr. Michael D. White 


Develop High-Performance Time-Domain Computational Electromagnetics Capability For 
RCS Prediction, Wave Propagation In Dispersive Media, and Dual-Use Applications 

The principle objectives of NCC3-885 

• Enhancement of an electromagnetics code (CHARGE) to be able to effectively model 
antenna problems. 

• Utilize lessons learned in high-order/spectral solution of swirling 3D jets to apply to 
solving electromagnetics project. 

• Transition a high-order fluids code, FDL3DI, to be able to solve Maxwell’s Equations 
using compact-differencing. 

• Develop and demonstrate improved radiation absorbing boundary conditions for 
high-order CEM. 

• Extend high-order CEM solver to address variable material properties. 

Enhancement of Electromagnetics Code CHARGE 

Charge is a structured finite-volume electromangetic solver that was part of an Electromangetic 
Grand Challenge Project. The solver originally utilized a 2nd-ordcr explicit Runge-Kutta time 
integrator and 3rd-order Van Leer flux vector splitting. One of the enhancements made to the 
code was improvement of the reflecting boundary condition at a conducting surface boundary. 
The improved boundary condition, while being more accurate, exposed a stability stiffness of the 
Runge-Kutta solver as originally implemented. The original solver was only stable due to the 
dissipation in the 3rd-order Van Leer, which was unsuitably small for typical problems. To 
mitigate this, the 2nd-order solver was replaced by a 4th-order Runge-Kutta, which is known to 
have much greater stability for hyperbolic problems. While the 4th-order Runge-Kutta has twice 
as many operations as a typical 2nd-order solver, it has more than twice the stability bound of the 
typical 2nd-order solver, making it more efficient. A comprehensive study of the efficiency of 
the parallel implementation of Charge was undertaken and published in collaboration with Dr. 
Jose Camberos of AFRL in the Center of Excellence for Computational Studies [1], Figure 1 
shows a generic finned missile being illuminated by a electromagnetic wave using Charge, while 
Figure 2 shows the scalability of the code. 

Another enhancement made to Charge was to add the ability to model patch antennas by use of 
equivalent currents [3, 4], This methodology allows for the patch characteristics to be determined 
separately, either by experiment or computationally, from the calculation on the body itself. 
Figure 3 shows a typical micropatch and its corresponding equivalent fields, while Figure 4 
shows radition pattern from the patch on a thin conducting plate. 
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(a) Antenna characteristics. 


(b) Graphical representation of the TM (0 .i) mode 
[Ey field]. 


Figure 3. Micropatch antenna. 


Previous work on a hybrid high-order/spectral scheme for analyzing round swirling jets was 
completed in early 2001. The formulation was a low Mach number approximation of the Navier- 
Stokes equations which is equivalent to a viscous formulation of the equations Lighthill used to 
analyze sound excitation from turbulence [4], This worked addressed topics of general interest 
for high-order methodologies in general, including various aspects of filtering for stability of the 
solution and efficient methods to parallelize these schemes [4, 5], The code utilized 6th-order 
compact finite-differences in the radial and axial directions hybridized with a (pseudo- )spectral 
method for the spectral direction. Temporal integration was achieved with either a 3rd-order 
explicit or a modified 2nd/3rd-order semi-implicit Runge-Kutta. An example of the cone- 
breakdown of a swirling jet is shown in Figure 5. 
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Angle (degrees) 

Figure 4. Comparison of radiated field with an integral method for 3.3 GHz patch on flat plate ( 1 8” x 1 8” x V4”). 




(a) Constant vorticity magnitude iso-surfaces at 
22 . 5 % maximum. 


(b) Cross-section of vorticity magnitude. 


Figure 5. Cone breakdown of a swirling jet: Re=880. S= 1 .33. 


As compact-difference schemes are notoriously unstable, explicit 8th-order filters were used to 
guarantee stability. For the spectral solution, it was discovered that the stability was strongly 
dependent upon the mode number divided by the radius as one approached the axis. Trying to 
solve the spectral terms implicitly did not improve stability as the angular and radial terms are 
actually dependent on one another [4, 5], To prevent the need for an extremely small time-step, a 
shaip filter was applied to the 0-direction that eliminated all wavelengths shorter than 2Ar, where 
Ar was the spacing in the radial direction at the axis. 
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Because the compact-differences are solved implicitly, parallelization is not, by nature, efficient. 
To overcome this difficulty, the compact-differences were broken up by either high-order 
explicit or one-sided compact-differences. This break up of the differences allowed for 
computational efficiency to remain above 50% for up to 64 processors. Similar runs without the 
breaking of the serial dependence had an efficiency of under 10% at 64 processors, as is shown 
in Figure 6. 



Figure 6. Efficiency of parallelization. * and * are SP3 compact-difference runs with serial dependence intact for 
grid sizes 101x131x7 and 101x131x17 respectively. ® are SP3 runs with explicit differences everywhere 
(101x131x7), while ■ and A are SP3 compact-difference runs with serial dependence broken by compact 
differences (for 101x131x7 and 101x131x17. respectively). • is the 101x131x17 compact-difference runs with 
serial dependence broken by compact differences run on the SGI Origin 2000. 


The knowledge from the high-order fluids code was instrumental in the transition of the high- 
order government code FDL3DI from a fluid solver to solving Maxwell’s Equations [6], These 
modifications are in continuing development under the working name AFRLCEM. This code 
currently has the ability to be run with up to 6th-order compact-differences spatially and utilizes 
a 4th-order Runge-Kutta for temporal integration. Implicit filters of up to lOth-order may be used 
for stability with minimal impact on the fidelity of the solution. Based on work by Drs. Miguel 
Visbal and Datta Gaitonde, both of AFRL’s Center of Excellence for Computational Studies, 
new absorbing boundary conditions were applied to Maxwell's Equations for the far-field 
radiation condition. Previously, high fidelity boundary conditions for Maxwell’s Equations 
required use of the technique known as Perfectly Matched Layers (PML). However, PML has the 
disadvantage of being difficult to implement and has known stability problems. The alternative 
boundary condition is tentatively known as the ETA boundary condition, where ETA is an 
acronym for Energy Transfer and Annihilation. This boundary is simple to implement and 
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robust. The idea behind it is to simply introduce a strong stretching at the absorbing boundary 
and lower the order of the filter as the boundary is approached [7]. The effect of the strong 
stretching is to reflect (transfer) part of the outgoing energy into high-frequency waves which are 
then annihilated by the high-order filters [7], As the filters are already being used for stability, 
there is negligible impact on solution time and no difficult coding is necessary. Figures 7 and 8 
show some of the results of the technique. 


Absortms rscon 




(a) No filtering. (b) Using implicit filter 

Figure 7. Propagation of a Gaussian Pulse from constant spacing onto a geometrically stretched mesh. 



(a) Computation (top half) compared 
to analytical solution (bottom half). 
Black line delineates absorbing 
boundary region. 


(b) Grid showing stretched (c) Solution comparison along the line 
absorbing region (in red) for at 0=45°. 

scattering problem. 


Figure 8. Scattered Ez field for a PMC cylinder at ka=50.67. 


In the past year emphasis was placed on the modeling of materia! properties, especially in high- 
order schemes [8]. One of the issues facing high-order schemes is the jump in material properties 
at a discontinuous interface. In generating the stencils for high-order schemes, it is generally 
assumed that the solution is continuous in the same order of spatial derivatives as the scheme. 
However, at a material interface, only the tangential E and H-fields and the normal D and B- 
fields are continuous. Likewise, the spatial derivatives of these fields are discontinuous at a sharp 
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interface. One advantage in treating these fields is that the jump conditions at the interface can be 
calculated for all derivatives. Because these fields have to be continuous for all time, the jump 
conditions in the derivatives can be calculated by matching Maxwell’s Equations on either side 
of the material interface. This work was the first to incorporate this information into the implicit 
compact-difference framework. 

Figure 9 shows a typical discontinuous derivative for a Gaussian pulse passing through a jump in 
material properties. Figure 10 shows the convergence of the derivative and filter in 1-D. As can 
be seen, treating the interface to 3rd-order was sufficient to retain 4th-order convergence 
globally. While technically, the filter should also be modified for the jump, the error introduced 
by filtering is found to be much smaller than that introduced by solving for the derivative; thus 
for practical purposes, it was found that no qualitative change in the solution was seen by using 
the original filters with the modified derivatives. In multiple dimensions, the solution is a bit 
more complicated, as the tangential derivatives must be calculated to properly treat the jump in 
the normal direction to the interface. Figure 1 1 shows the effect of modifying the 1-D solution by 
using the tangential derivatives for a single calculation of the time derivative of the E-field. 
Figure 12 shows the comparison of a material coated perfect magnetic conducting (PMC) 
cylinder illuminated by a plane wave. 



Figure 9. Computed discontinuous derivatives at a material interface. 
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Derivative treatment at interface 


Filter treatment at interface 




Ln(Ax) 


(a) Convergence of derivative. (b) Convergence of filter. 

Figure 10. Gaussian pulse passing through a material interface with e r = 4 and |p = 2. 



r 


Figure 11. Effect of solution with and without modification at interface for a single calculation of right hand side for 
a material coated PMC cylinder. Material properties between PMC cylinder and free-space (0.5 < r < 0.821 ) are e/e 0 
= 2.2, p/po=4. 1,0 = 45°. 
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Figure 12. PMC cylinder with coating: g r = 1, e r = 4.1 illuminated by a TM pulse. The cylinder has radius V 2 , and 
the material layer has thickness r m = 0.3124. The period of the incident pulse is t = 0.525 and the total field zone 
(TFZ)/scattered field zone (SFZ) is located at approximately r = 1.3. The top half of the figure is the computed 
E z -field solution, while the bottom half is the series solution. 

In collaboration with Dr. Brad Edgar of Arkansas State University, an effort was recently begun 
to look at simulating dispersive media with these high-order schemes. Dr. White was 
instrumental in advising Dr. Edgar on the simulation of Debye and Lorentz media. The media 
was coupled to Maxwell's Equations by use of auxiliary differential equations to solve the 
frequency dependent response. Dr. White has subsequently incorporated a slightly different 
auxiliary difference formulation into the AFRLCEM framework. Initial results of this latter work 
can be seen in Figure 13 which shows a Gaussian-modulated pulse traveling through a tenuous 
plasma. This work on dispersive media is continuing under a different grant. 




(a) Initial condition. (b) Solution at t= 1 4.5 . 

Figure 13. Gaussian-modulated pulse propagating through a plasma. At =0.005 and the non-dimensional plasma 
parameter. ayi = !0n. Initial pulse splits into two pulses moving in opposite directions. 
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Mr. Gerald Trummer 


Mr. Gerald Trummer is Systems Analyst for the Computational Sciences Branch and the 
Computational Sciences Center of Excellence at the Air Force Research Laboratory’s 
Aeronautical Sciences Division at Wright-Patterson AFB. He is responsible for integrating the 
hardware and software of the Center’s heterogeneous computer environment into a homogeneous 
resource for the high-speed calculation and visualization of computational physics data. He 
analyzes the computational needs of the Scientists and Engineers of the center and pursues 
meeting those needs in the acquisition and development of hardware, software and visualization 
tools. He supports the S&Es in the use of computer graphics software and hardware, scientific 
visualization tools, and parallel computation. He maintains and administers the Aeronautical 
Science Division’s cluster of high-end graphics workstations and the visualization center. In 
addition, he administers a 32-node and a 13-node cluster of linux workstations for parallel and 
distributed computation of numerical problems. He is also responsible for the production of 
scientific visualization movies representing the work of the S&Es for presentation at conferences 
and lectures. Mr. Trummer also acts as scientific visualization liaison between the Air Vehicles 
Directorate and the Aeronautical Systems Center’s High Performance Computing Center and in 
this capacity, works with the ASC HPC’s visualization experts to obtain and improve resources 
and services available to the S&Es of the Directorate. 

Mr. Trummer installed hardware and software upgrades on existing AFRL/VAA UNIX 
workstations, which included operating system upgrades, security patches and graphics 
applications. Integrated new high performance workstations into the existing AFRL/VAA 
workstation network cluster, which included installing hardware, operating system upgrades, 
security patches, graphics applications, configuring system files and adding user accounts and 
files. Installed QUALSTAR TLS 4210 TAPE LIBRARY and software to perform 
comprehensive tape backup of all user data on the AFRL/VAA workstation cluster. Integrated 
upgraded workstation and video hardware and software in the AFRL/V AAC conference room to 
perform scientific data acquisition from wind tunnel test results. Made scientific visualization 
movies of CFD and CEM data for the scientists of the Aeronautical Sciences Division for 
presentation at AIAA conferences, symposiums and lectures. Among which included Dr. 
Gordnier’s “Computation of Three-Dimensional Non-linear Panel Flutter” and Dr. Rizzetta’s 
“Supersonic Compression Ramp Flow” for presentation at AIAA conferences, symposiums 
and/or lectures. 

Mr. Trummer upgraded the capability and capacity of the comprehensive data backup of the 
AFRL/VAA workstation cluster. Installed a new SGI Origin 300 server for high performance 
CFD and CEM calculations including hardware for mass file storage and software for post 
processing visualization of scientific data. Installed new hardware and software in support of the 
Computational Sciences Center of Excellence Immersion room, including an Accom Digital 
Disk Recorder for the creation of CFD and CEM visualization movies for presentation in the 
Immersion room and remote sites. 

Mr. Trummer integrated high performance graphics workstations into the existing network 
cluster which included installing new hardware, configuring operating system software. 
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installing patches and graphics applications, adding user accounts and files, and testing and 
evaluating the security of the installation. 

Mr. Trummer received a Notable Achievement Award during the period April 1 to June 30, 2001 
for providing exceptional computer support and services to the Computational Sciences Branch. 
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